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ABSTRACT 

We analyse the effect of two relevant physical constraints on the mass multiplicity 
function of dark matter halos in a Press-Schechter type algorithm. Considering the 
random-walk of linear Gaussian density fluctuations as a function of the smoothing 
scale, we simultaneously i) account for mass semi-positivity and ii) avoid the cloud-in- 
cloud problem. It is shown that the former constraint implies a severe cutoff of low-mass 
objects, balanced by an increase on larger mass scales. The analysis is performed both 
for scale-free power-spectra and for the standard cold dark matter model. Our approach 
shows that the well-known "infrared" divergence of the standard Press-Schechter mass 
function is caused by unphysical, negative mass events which inevitably occur in a 
Gaussian distribution of density fluctuations. 
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1 INTRODUCTION 

The problem of understanding the origin and the evolution of the density fluctuation field represents one of the fundamental 
issues of modern cosmology. Primordial perturbations are believed to be generated from vacuum fluctuations of a weakly 
coupled scalar field, during an inflationary stage. Because of this reason, it is often assumed that the cosmological density 
fluctuations, during their linear evolution, made up a Gaussian random field. This assumption might be also motivated on the 
basis of the Central Limit Theorem. It is well known, however, that this choice violates the mass semi-positivity requirement 
(5 > — 1). In fact, a Gaussian field with zero mean always admits a finite probability of assigning events with values of the 
random variable lower than —1; furthermore, this probability raises when the field variance increases. Thus, if one applies a 
coarse graining procedure to a hierarchical Gaussian density field (defined by da 2 /dRf < and a 2 — > +oo, when the filtering 
length Rf — + 0), one notes that the probability of finding regions with "negative mass" grows with decreasing Rf. Therefore, 
we expect this feature of the Gaussian distribution to affect the counting and arrangement properties of low-mass objects. 

In this paper we study the effects caused by avoiding this shortcoming of the Gaussian choice in obtaining the dark 
halo mass function for a hierarchical scheme of structure formation. Note that allowing for the physical constraint deriving 
from mass semi-positivity amounts to introducing a sort of minimal non-Gaussianity in the density fluctuation field. Some 
generic features of the mass function obtained by non-Gaussian perturbations have been already investigated by Lucchin & 
Matarrese (1988), who found a general tendency for an increase of the number of high-mass objects, and by Sheth (1995). A 
related problem has been studied by Catelan et al. (1994), who analysed the effects of forbidding negative mass events on the 
two-point correlation function of regions up-crossing a density threshold. 

We know that hierarchical clustering proceeds from the "bottom-up": low-mass halos form first while bigger ones are 
created by aggregation and merging of still existing objects. For this reason, those models which aim at determining the mass 
function of cosmic structures in a hierarchical aggregation scenario must afford the problem related to the existence of sub- 
condensations inside clumps on larger scale. This "cloud-in-cloud" problem represents the main drawback of the classical 
Press-Schechter theory (Press & Schechter 1974) in which, on the grounds of the spherical collapse model, the virialized 
objects are identified with those regions where the filtered density field, during its linear evolution, becomes greater than a 
critical threshold, corresponding to a certain density contrast S c - 



2 C. Porciani, F. Ferrini, F. Lucchin and S. Matarrese 



Our approach is an extension of the one developed by a number of authors (Peacock & Heavens 1990; Cole 1991; Bond 
et al. 1991), which succeeds in excluding sub-condensations from the count of bound objects. [Alternative approaches to the 
mass function have been followed by Manrique & Salvador-Sole (1995) and Cavaliere et al. (1995).] According to this method, 
in each point one considers the "trajectory" 8(Rf) of the density field as a function of the filtering radius and determines 
the largest Rf (and then the largest possible mass) at which 8(Rf) crosses the threshold S c . This is enough to solve the 
cloud-in-cloud problem: all the objects selected in this way cannot have been included in bigger condensations, since the 
surrounding regions, filtered on all the larger scales have density smaller than the critical one. Moreover, the objects which 
could be associated to threshold crossings occurring on smaller scales must not be considered, since the collapse of a structure 
erases every sub-condensation. In much the same way we succeed in avoiding the contribution from negative mass events 
to the count of dark matter halos by imposing a boundary at <5„ = —1 to the random walk S(Rf). This constraint has two 
important and complementary effects: it implies a substantial decrease of the number of low-mass objects, thereby eliminating 
the low-mass divergence of the Press-Schechter mass function, which is balanced by an increase of the number of high-mass 
objects. Moreover, also the redshift evolution of the resulting mass function is largely modified compared to the standard 
Press-Schechter theory. 

We mostly follow the formalism by Bond et al. (1991); in Section 2, we derive the mass function according to their 
prescription while, in Section 3 (but see also the appendix), we modify it to account for the simple but physically relevant 
non-Gaussian feature discussed above, which will be shown to have a strong impact on the low-mass behaviour of the halo 
multiplicity function. Section 4 contains a brief discussion and some conclusions. 



2 THE MASS FUNCTION OF DARK MATTER HALOS: A STOCHASTIC APPROACH 

In this section we briefly review the mathematical formulation of the "random-walk" approach sketched above. After some 
general remarks, we focus on the configuration that considers a sharp fe-space filter in the presence of a single absorbing 
barrier set at the threshold density. We basically follow the approach by Bond et al. (1991), although some aspects of the 
formulation are slightly modified. 

Let us assume that the primordial density fluctuations form a homogeneous and isotropic Gaussian random field <5(x, z), 
uniquely specified by its power-spectrum P(k, z) (in the following, the power-spectrum at z = will be simply denoted by 
P(k)). If one identifies the collapsed regions with those points where the filtered mass density field lies above a constant 
threshold, one can allow for the redshift dependence of 8 in terms of the growing mode of linear perturbations, D + (z). 
Following Bond et al. (1991), however, we ascribe the redshift dependence to the formation threshold and consider the density 
fluctuations as a static random field, <5(x), normalized to its linear extrapolation to the present time. The evolution of 8 C is then 
fixed by the background cosmology: the spherical collapse model gives 8 c (z) — A(z)/D + (z) where the function A(z) depends 
weakly on the values assumed at redshift z by the density parameter Q, the cosmological constant A vac and the Hubble 
constant H (Lilje 1992). In a matter dominated Einstein-de Sitter universe (the only case considered in this paper) the 
threshold increases with redshift according to the relation 8 c (z) — 8 c /D + (z), with D+(z) = 1/(1 + z) and 8 C — const = 1.686. 

By convolving the density field with the filter function VK(|x — x'| , Rf) one obtains a new field <5(x, Rf), which is defined 
in the four-dimensional space (x, Rf) and which, in general, is not homogeneous and isotropic along the Rf direction. The 
dependence of 5(x, Rf) on Rf can be deduced with complete generality using the Fourier transform of the density field, 
5(x, Rf) = (2-7r)~ 3 J S(k)W(kRf)e~'^ i x d 3 fc. In fact, an infinitesimal change of Rf affects the value of <5(x, Rf) according to 
the relation 

W(x,«/)_ 1 fsj^/fK-^k^^Rf). (1) 



dR f (2tt)3 J w dR f 

Due to the stochastic nature of <5(x) and to the linearity of Eq. (1), it follows that r?(x, Rf) is also a zero mean Gaussian 
random field, which is therefore uniquely defined through its auto-correlation function 
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obtained by using the definition of power-spectrum (<5(ki)<5(k2)} = (27r) 3 5i)(ki + k2)P(fci), where fo(y) represents the Dirac 
delta function, and by defining r = |xi — x 2 |. 

The equality that defines r\ has the form of a Langevin equation for the smoothed density fluctuation field changing under 
the action of the stochastic force 7j(x, Rf). Unfortunately, for the most popular filter functions, such as Gaussian and spherical 
top-hat, ?](x.,Rf) becomes a "coloured" noise, whose correlation properties along the Rf direction make the problem too 
involved to be afforded by analytical means. It has been shown by Bond et al. (1991) that the problem becomes much more 
tractable if one smooths the density field by a "sharp fc-space" filter, i.e. W(k, kf) = 9{kf — k), where 8(x) is the Heaviside 
step-function and kf oc 1/Rf- The calculation of the mean mass enclosed in the filtering volume has been performed by 
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Lacey and Cole (1993) and gives M(kf) = 6ir 2 {g)/k 3 f. With such a filter, decreasing the radius corresponds to adding up a 
new set of Fourier modes of the unsmoothed distribution to S(Rf); for a Gaussian random field, this increment is completely 
independent of the previous steps, so that the trajectory S(Rf) represents a Brownian motion. 
In practice, one can use the variable kf as "time" variable, obtaining 



cW(x, kf) 



(2tt)3 / 



S(k)S D (k f - k)e 



(3) 
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where £(x, kf) is a new Gaussian stochastic field. By averaging over the statistical ensemble one then finds (£(x, kf)} = and 



(C(x 1 ,fc /1 )C(x 2 ,fc /2 )) = -l ¥ fc^ 1 P(fc /1 ) Sin(fe/ir) fo(fc /1 
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kf 2 ). 



In an arbitrary point of space the density fluctuation field evolves with kf according to the Langevin equation 

ds(kf) 



dkf 



= £(*/), 



(4) 



(5) 



with the initial condition 5(0) = 0. By averaging 8(kf) = f^£,{s) ds over the statistical ensemble, one obtains (S(kf)) = and 
(5(kfi)5(kf2)) = (1/2-k 2 ) J™ 11 '*/ 1 ''/ 2 ' s 2 P(s) ds, which completely determine the probability density W(5, kf) of the Gaussian 
process S(kf), namely a zero mean Gaussian distribution with variance a 2 {k f ) = (1/2tt 2 ) J Q f s 2 P(s)ds. These results show 
that our physical system is dynamically equivalent to a set of particles undergoing one-dimensional Brownian motion x(t) 
with diffusion coefficient varying with time. This analogy becomes even more evident if one identifies the time variable with 
the variance A = u 2 (kf) of the filtered density field. In such a case, the stochastic process looses the time-dependence of the 
diffusion coefficient and becomes a Wiener process, 85(A)/dA = £(A), with (C(A)) = and (£(Ai)£(A2)) = <5d(Ai — A 2 ). 

So far we have analysed the problem described by the Langevin equation (5) plus "natural" boundary conditions: 
lims^ioo W(5, kf) = 0. Our aim is, however, to compute the fraction of all trajectories which have crossed, at least once, the 
threshold for structure formation at a given resolution kf. Such a quantity can be evaluated by putting an absorbing barrier 
in 5 = 5 c (z): when a realization of the random process S(kf) reaches for the first time the level 5 C one stops to count it, so 
that one always knows how many realizations have not yet reached the barrier. In order to analytically solve this problem it 
is more convenient to follow directly the behaviour of the probability density W(5, kf), which can be easily shown to obey the 
Fokker-Planck equation, 



dW(5,A) 



ld 2 W(5, A) 
2 



(6) 



OA 2 dS 2 

The solution of equation (6) with absorbing boundary condition in 5 = 5 C , W(<5 C ,A) = 0, and with the initial condition 
W(5, 0) = 5d(8) has been obtained for the first time by Chandrasekhar (1943); it reads 

S 2 \ ( [5 - 25 c 



W(S, A, 5 C 



exp 



2A 



exp 



2A 



(7) 



By integrating the previous expression over the allowed region, one obtains the probability that, by the "time" A, a 
particle has not yet met the barrier, S(A, 5 C ) = J_° W(S, A, 8 C ) dS. Then the probability that, during its stochastic motion, a 
given trajectory has crossed the critical level at a variance lower than A can be deduced from the probability conservation law, 
Q(A, 5 C ) = 1 — S(A, 5 C ). By differentiating with respect to A one obtains the probability distribution function of first-crossing 
variances, 



/(AA) = 



dQ(A,S c ) 
dA 



d_ 

dA 



r 
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W(<5, A, <5 C ) d5 = 



18W(5,A,S C ) 
~ 2 
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Bond et al. (1991) used these results to get an improved, Press-Schechter-like expression for the mass function free of 
the cloud-in-cloud problem. The function /(A, 5 C ) dA yields the probability that a realization of the random walk is absorbed 
by the barrier during the time interval (A, A + dA), or, thanks to the ergodic theorem, the probability that a point is involved 
in the collapse of a structure in the mass range [M(A + dA), M(A)]. The comoving number density of structures with mass in 
the range (M, M + dM) present at redshift z is therefore given by 



n(M,z)dM=^f(A,5 c (z)) 



dA 



dM 



dM, 



which, using equation (8), becomes 
n{M,z)dM= { - Q)5c[l + z) 1 



2tt 



M 2 AV2(M) 



dlnA 



dlnM 



exp 



5 2 (l + z) : 
2A(M) 



dM. 



(9) 



(10) 



This equation is identical to the Press-Schechter formula, including the well-known "fudge factor" of two. 
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3 THE ZERO MASS BARRIER 

The technique adopted to solve the cloud-in-cloud problem hides an inconsistency that is peculiar to every model which tries 
to describe the density fluctuations by means of a Gaussian field. Indeed, to deal with a universe which does not contain 
regions of negative mass, one must assume that the density field £>(x) is semi-positive definite and, as a consequence, that the 
corresponding fluctuations <S(x) get only values larger than or equal to —1 everywhere in space. Moreover, if one considers 
a window function that is semi-positive definite (except for a set of measure zero), also the filtered density fluctuation field 
must only get values 5(Rf) > 5 V = — 1. All this conflicts with the Gaussian hypothesis which in any point predicts 8(Rf) < — 1 
with finite probability. 

We want now to modify the algorithm presented in the previous section to account for this physical constraint. Unfortu- 
nately one cannot easily forsake the Gaussian assumption, because of its relevant role in the construction of the Fokker-Planck 
equation (Pawula theorem; see, e.g., Risken 1989), so one should devise a stratagem able to provide a theory equivalent to a 
non-Gaussian one, but which remains tractable. 

The existence of regions with negative mass (8 < —1) in the Press-Schechter approach is caused by two main reasons: i) 
the initial conditions are represented by a Gaussian field; ii) the perturbations are evolved within the linear approximation. 
While the first motivation represents a real inconsistency if applied to the density fluctuation field, the second one is much 
more subtle: we know that the real evolution of the perturbation field differs sensibly from its linear approximation. Following 
the exact dynamics and starting from a self-consistent distribution of the density field (i.e. one not containing negative 
mass events) one would never obtain regions with negative mass. However, in order to make quantitative predictions, one 
extrapolates the validity of the linear theory beyond its limits and, keeping in mind the spherical model, one assigns a rule- 
of-thumb for modelling the collapse of overdense fluctuations. Similarly, one should remind that underdense fluctuations 
also experience non-linear evolution, so that one could devise a suitable short-cut to account for the non-linear dynamics 
of underdense regions. In this section we treat this problem in a formal way: we explicitly forbid our random trajectories 
to enter the unphysical region, by putting a reflecting barrier in 8 — 8 V = — 1. To account for the time dependence of the 
density fluctuation field, we assume that the value 8 V increases with redshift according to 5 v (z) = 8 v /D + (z). The aim of this 
approach is to understand which features of the mass-function actually depend on the presence of negative mass events in 
the probability distribution. 

The possibility to obtain an analytical solution is once again restricted to the choice of the sharp fc-space filter, which 
is, however, slightly inconsistent with the semi-positivity assumption for the filter function. In fact, the sharp cut in Fourier 
space implies, by the indeterminacy principle, an infinite series of oscillations in configuration space where the filter assumes 
many times negative values. On the other hand, the absolute value of this function decreases quite rapidly as one goes away 
from the origin, so that the integral which defines the filtered density field is largely dominated by regions where the filter 
is positive. We may then expect that the error one makes by using this window function is small (notice, moreover, that a 
similar problem is present for the absorbing barrier set at 8 C ). 

The choice of a sharp fc-space filter together with the presence of a reflecting barrier in 8 — S v and of an absorbing one 
in 8 = S c allows to write the Fokker-Planck equation (6) for the probability density W(S, A), with the boundary conditions 



2 88 



= 0, W(<5 C ,A) = (11) 



and the initial condition W{8, 0) = 5d(8). Here J (8, A) represents the probability density current and the equation which 
involves it characterizes the reflecting barrier. 

In this case the evolution of the system is analogous to that of a set of particles undergoing Brownian motion in the 
presence of the potential V(x) = +oo if x < 5 V , V(x) = const if 8 V < x < S c and V(x) — — oo if x > S c . Starting from this 
analogy, one may wonder for which properties of the first-crossing distribution one should expect relevant differences from 
the case previously considered. Indeed, it is easy to deduce that the first-crossing times are smaller on the average, since the 
reflecting barrier reverses the motion of those particles which hit it, forbidding their dispersion to very large distances away 
from the absorbing boundary. As in our analogy the time variable corresponds to the variance of the filtered density field, one 
should expect that the reflecting barrier increases the number of crossings at small variance; in practice, the numerical density 
of small mass objects should decrease while that of large mass clumps should increase. On the contrary, those realizations 
that reach the critical level in a very short time, describing a "quasi-coherent" trajectory headed forward, are not influenced 
by the presence of the reflecting barrier. One should then expect the numerical density of very large mass objects to be 
unaffected by our procedure. A numerical computation of the mass function performed by using the adhesion approximation 
in one dimension seems to share the same behaviour (Williams et al. 1991). However, we stress that N-body simulations of 
hierarchical clustering in three dimensions show better agreement with Press-Schechter results (Lacey & Cole 1994). 

To solve the Fokker-Planck equation let us call x the random field 8 and t the independent variable A. Call then A c the 
position of the absorbing barrier and —R v that of the reflecting one. As initial condition one assumes a Dirac delta function 
set in x = 0. To simplify the equations one can shift the origin so that it corresponds to the location of the reflecting boundary; 
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»(M,,) = ^ 

2(l + <5 c ) 2 (l + z) 2 M 2 



din A 



dlnM 



in such a way the absorbing barrier is found in x — A c + R v and all Brownian particles leave from x = R v . One has then to 
solve the equation 

dW(x,t) _ 1 8 2 W(x,t) 

dt 2 dx 2 1 ' 

together with the boundary conditions dW(0,t)/dx = 0, W(A C + R v ,t) — and the initial condition W(x,0) = 5d{x — R v ). 
The problem can be solved by separation of variables. One finds W(x, t) = 2 J^^Lo 4>n(x)4>n(Rv) exp(— |A n i), with eigenvalues 
An = [(2n + l)n/2(A c + Rv)] 2 and orthonormal eigenfunctions <f>„(x) = (A c + Rv)^ 1 ^ 2 cos (\/X^x). The probability density 
function then reads 

MM, - ^ f - (&±£>«.) ~ (^«) «P (-g^.) . 

while the first-crossing rate across the absorbing barrier in x = A c is 

n— x 7 x 7 

To obtain the mass function deriving from the above crossing rate one must go back to the original physical variables S 
and A and replace /(A, 5 c {z)) = T C (A) in Eq.(9). Assuming fi = 1 one gets 

v\ „ /(2n + l)7r\ / (2n + l) 2 n 2 A(M)\ 

( + } IW) J ° XP (- 8(1 + ^(1 + ^ J ■ < 15 > 

We can now compare our result with the Press-Schechter mass function for scale-free power-spectra P(k) = Ak np , 
where n v is the primordial spectral index (let us remind that detailed analyses of the mass function of clumps in A^-body 
simulations with scale-free initial conditions have been performed by Efstathiou et al. 1988 and by Lacey and Cole 1994). In 
such a case, selecting a sharp fc-space filter, one has A(M) = a 2 (M,z = 0) = (M/M )~ 2a with M = &-k 2 {q) (A/12-ir 2 a) 1/2a 
and a = (n p + 3)/6. Replacing these expressions in Eq.(15), we obtain 

■>-¥#(=&)) " i?- 1 "- + - (s) - (-<* - (to) . 

where Mr(z) = Mo [7r 2 /8(l + S c ) 2 (l + z) 2 ] 1 ^ 2a represents the characteristic mass of the distribution at redshift z. Since the 
function M 2 n(M, z) depends only on the ratio M/Mr(z), our solution keeps the self-similarity property characterizing the 
Press-Schechter theory. 

In Figure 1 we show the behaviour of our new mass function by plotting M%(z)n(M, z) / (g) which depends only on the 
ratio M/Mr and is not affected by the normalization of the power-spectrum. We consider different values of the spectral 
index (n p = —2, — 1, 0, 1) comparing our solution with the Press-Schechter one. It is evident that these mass functions are 
very different in the low-mass tail while the discrepancy tends to disappear for large masses. In fact n(M) defined in Eq.(10) 
diverges as M — > whereas ours goes to zero in the same limit. 

This behaviour is quite interesting: in fact it is rather unnatural to imagine that n(M) grows unbounded for M — > 
in a hierarchical scenario where at any time aggregation processes are able to conglomerate small objects. Indeed, a study 
of the time evolution of our mass function (Figure 2) shows how hierarchical clustering displaces power from small to large 
scales. The presence of a peak in the mass function confirms the existence of a time-dependent characteristic mass. This peak, 
however, becomes less and less prominent as time goes on. It is easy to show that the maximum value of the mass function 
is reached for M/M R ~ (a/a + l) 1/2a = [{n p + 3)/(n p + 9)] 3/(n " +3) . 

We can now apply Eq.(15) to a physically sensible model, such as the standard Cold Dark Matter (CDM) scenario. We 
compute the appropriate A(M) by using the transfer function given by Bardeen et al. (1986), with tlx = 1 and Ho = 50 km 
s _1 Mpc -1 , and by normalizing the spectrum so that the mass variance is unity in a top-hat sphere of radius 8 hT 1 Mpc. 
As Figure 3 shows, also in this case the mass function presents a strong low-mass cutoff. In fact, the number density of halos 
reaches its maximum at M ~ 4.67 x 10 11 Mq and is strongly damped at much smaller masses. We will examine the physical 
reliability of this cut-off in the following section. 

To understand how all the mass is divided among the various objects, in Figure 3 we plot also the "multiplicity function" 
M 2 n(M)/(g) which gives the mass fraction contained by halos in unit range of InM. Obviously, the dlnM integral of this 
dimensionless distribution, performed over the whole mass spectrum, gives unity. 

In Figure 4 we show the time evolution of our mass and multiplicity functions; due to the characteristic scales inherent in 
the CDM power-spectrum the self-similarity property is now lost, even though the mass once again flows from small objects 
to bigger ones. In order to quantitatively follow this process, we study the time dependence of the typical cluster mass of the 
distribution as defined in the kinetic theory of aggregation (see, e.g., van Dongen & Ernst 1988) 
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»„ ^ f °° M2n (M,z)dM /■« M 2 n(M,z) 

M ci («) ee ^ = / ,\ dM. (17) 

J™Mn(M,z)dM J (g) 

In the interval < z < 2, our results are approximately described by the power-law 

M cl (z) = M ci (0)(l + Z )-' 3 (18) 

with M c ;(0) ~ 6.39 x 10 14 M Q and /3 ~ 3.36; however, a deeper analysis reveals that the quantity (3 e ff(z) = I ^'h^i+i^ I s l° w ly 
increases with z. As one would expect, we find (3 e ff(z) ~ l/a e ff(z) where a e ff(z) = — ''^y |m c; (z)- We stress that for 
the Press-Schechter formula this equality does not hold because the integral that defines M ct (z) comes out mostly sensitive 
to low-mass abundances (Colafrancesco, Lucchin & Matarrese 1988). 

According to high resolution iV-body simulations, for 10 10 Mq < M < 10 13 Mq the mass function of a standard CDM 
scenario is well described by a power-law of index —2.0 ± 0.1 in a wide redshift range (Brainerd & Villumsen 1992). This 
feature is not displayed by the Press-Schechter solution whose logarithmic derivative assumes values ~ —1.8 in the same 
interval. So, in order to test our solution, we study the functional dependence of its slope on mass. We find that only for 
z >2 our solution behaves like M~ 2 , with very good approximation in the mass interval under consideration. Anyway, we 
notice that the choice of the values of the parameters 5 C and v — M (k /)fc 3 / '(g) plays a fundamental role in this kind of 
comparison. We remind, for example, that the value S c = 1.686, obtained from the spherical collapse model, should not be 
preferred when one uses a sharp fc-space filter. Probably, as suggested by Williams et al. (1991), the wisest choice is to use 
numerical simulations to select S c and v as best fitting parameters (see also Bond & Myers 1993, Ma & Bertschinger 1994, 
Klypin et al. 1995, Monaco 1995). It is known however that these best fit values depend both on the filter used and on the 
method chosen to identify the halos in the simulations (Lacey & Cole 1994, Gelb & Bertschinger 1994). 



4 DISCUSSION AND CONCLUSIONS 

In this paper we have derived a new analytical expression for the dark halos' mass function that develops in a hierarchical 
clustering scenario. This result has been achieved by simply modifying the Press-Schechter model to allow for mass semi- 
positivity. 

Technically speaking, in every point of configuration space we have studied the random walk of the coarse-grained density 
field as a function of the smoothing scale. Essentially, we have derived the probability distribution of the filtering lengths that 
characterize the events of first up-crossing of a critical threshold 8 C in the presence of a barrier set in 5 = 8 V = — 1. This 
boundary had to assure the reflection of the incident "random-walkers" to allow for mass semi-positivity. 

For power-law spectra we have obtained a mass distribution whose low-mass tail, for any z, has the general form 

n(M,z) ~ Ci(z)M- 2{a+1) exp(-C 2 (z)M- 2a ) M < M R (z), (19) 

where the characteristic mass Mr and the parameters Ci and C2 depend on the selected threshold and power index. Further- 
more, the mass function presents a maximum for M/Mr ~ [(n p + 3)/ (n p + 9)] 3 /( n p+ 3 ) and in the high-mass tail asymptotically 
reaches the Press-Schechter solution. During its time evolution the mass multiplicity function keeps self-similar. 

Working directly with many realizations of a Gaussian random field and performing an object by object analysis, a 
number of authors (Bond et al. 1991, Williams et al. 1991, White 1995) showed that for M « M,(z) = M /[5 C (1 + z)] 1,a the 
excursion set approach does not mimic with good approximation the kinetics of mass aggregation in a hierarchical scenario, 
while for M > M» the correspondence appears quite satisfactory. Since for many spectral indices, the mass for which our 
function reach its maximum value is of the same order of magnitude as M* , it would be interesting to test our solution against 
numerical simulations. 

A careful analysis of the mass function of scale-free models in an Einstein-de Sitter universe has been recently performed 
by Lacey & Cole (1994), using the high resolution P 3 M code of Efstathiou et al. (1985). Taking advantage of self-similar 
scaling, they succeeded in reducing Poisson fluctuations in counts by averaging the various outcomes obtained at different 
timesteps of the same simulation. In such a way they could amplify the available mass-range. Nevertheless, though they 
considered three different spectral indices (n = —2, — 1, 0), only for n = their low-mass limit is such as to make potentially 
observable the cut-off implied by our solution. Moreover, n — is the only index for which the different choice of the filtering 
method (top-hat with M(Rf) — 47r(£>)i? 3 c /3 in the work of Lacey and Cole, sharp-k with v = 6tt 2 in ours) has no effect. 
However, the numerical results agree quite well with the Press-Schechter solution showing no trace of a low-mass cut-off. For 
example, for n = and M/Mr ~ 0.2 (which is the lowest mass considered in the mass function obtained from the numerical 
data) our solution is almost a factor of 3 smaller than the Press-Schechter one which, on the contrary, overestimates the 
counts by ~ 30%. 

We are not surprised by this disagreement, since the method we used to allow for mass semi-positivity is a very crude one. 
We think that only by dealing with a field that has the correct statistical properties one can determine the exact position of 
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the cut-off. Our algorithm is only able to show that a low-mass cut-off must exist and to explain the origin of the divergence 
of the Press-Schechter solution for M — > 0, as being due to unphysical negative mass events in the Gaussian distribution. 

In order to test our solution also in the CDM case we refer to a recent paper by Efstathiou (1995), who analysed the 
redshift evolution (for 1 < z < 3 and M > 10 10 M o ) of the integrated mass function N(> M,z) = n(M',z)dM' in a 
highly biased, CDM dominated Einstein-de Sitter universe. By comparing our results with these data we find a fairly good 
agreement for z = 3, while we predict a larger number of counts (~ 75% more) for M ~ 10 10 Mq at z = 1. Actually, at high z 
we are testing the high-mass end of the distribution, where our solution is practically indistinguishable from the one of Press 
and Schechter, that fits the data quite well. On the contrary, at intermediate redshifts, the mass-range in analysis involves 
masses just above the peak, where our solution implies more objects than the Press-Schechter one (see also the figures 3 and 
4). 

Even though we faced the negative-mass problem in quite a formal way, our treatment allowed to understand the origin 
of the "infrared divergence" of the Press-Schechter theory. We do not claim that we are able to indicate the best solution, 
since we believe that the correct answer must rely on the intrinsic non-Gaussian nature of any density fluctuation field. In any 
case, we think that the different behaviour of our solution compared to the Press-Schechter one, if confirmed and improved 
by more detailed modeling, might be useful to develop a consistent picture for the formation of dark halos. In fact, it is 
known that, assuming a constant mass-to-light ratio, the Press-Schechter formula predicts too many low-mass objects with 
respect to the observed galaxy luminosity function (Bond et al. 1991). However, to deal with galaxies one should consider 
many astrophysical and hydrodynamical effects (gas cooling in non stationary conditions, star formation and so on) that are 
supposed to be fundamental issues of galaxy formation but are very hard to model. These subjects are clearly beyond the 
purposes of this work. 

In summary, we have shown that the low-mass divergence of the Press-Schechter mass function can be ascribed to the use 
of Gaussian fields to describe the cosmological density fluctuations, which assign a finite probability to events with a negative 
mass; since this probability comes out directly proportional to the variance, in a hierarchical clustering model the reliability 
of theoretical predictions should get worse as the variance increases i.e. as the mass decreases. We have shown that this is 
indeed the case. We believe that a reliable model able to make quantitative predictions should not leave truly non-Gaussian 
fields out of consideration. Even though there are good reasons to think that the primordial gravitational potential was very 
nearly Gaussian distributed, one should derive the statistical features of S from the non-linear fluid-dynamical equations. 
Only in this way one would be able to obtain the correct statistical properties of the density fluctuation field. We will return 
to this subject in a future work. 
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APPENDIX Al: VOID REGIONS AND THE MASS FUNCTION WITH TWO ABSORBING BARRIERS 

In Section 3, by modifying the excursion-set approach, we have seen how the introduction of a second barrier, able to limit 
the dispersion of the values assumed by <5(x, 7?/), would produce relevant changes in the behaviour of the mass function. We 
want here to investigate the effect of changing the nature of this new barrier, with the aim of simulating a particular physical 
situation, namely the existence of void regions. 

Let us consider a physical density fluctuation field, which is clearly larger than or equal to —1 everywhere in space; 
provided the applied filter is also semi-positive definite and correctly normalized, the only possibility to obtain the value — 1 
for the smoothed density fluctuation field is realized when the entire region which contributes to the convolution integral 
defining 5(x, Rf) is void. 

Let us then imagine to use a filter with radius Rf and find <5(x, Rf) — —1 in some arbitrary point x; for every smoothing 
radius smaller than Rf (therefore for every variance A > A(Rf)) the considered region must still be void. In practice one must 
obtain 5(x, A) = —1 for every A > A(Rf). Therefore, once the value —1 has been attained at the variance A(Rf), the density 
field in the point x cannot assume any other value: also the barrier set at 8 — 8 V behaves as an absorbing one. In such a way 
one succeeds in accounting for all the points included inside void regions. For each of these points x„, in fact, there exists a 
smoothing radius Rf corresponding to the minimum distance of the point to the boundary of the void region, such that, for 
any Rf < Rf, one measures 5{x v , Rf) = —1 (obviously, we are dealing with filter functions that do not vanish only in a finite 
region of space, as, for example, the top-hat one). 

Once again, in order to obtain quantitative results, we need to use the sharp fc-space filter even though this weakly violates 
our hypotheses. Then, one has to work out the Fokker-Planck equation (6) with the boundary conditions W(6 V ,A) — 0, 
W(5 C ,A) = and the initial condition W(8, 0) = 8d(8). Hence, using the same notation of section 3, one has to solve Eq.(12) 
with the boundary conditions W(0,t) — 0, W(A V + A c ,t) — and the initial condition W(x,0) = 8o(x — A v ) (the previous 
parameter R v has been replaced by A v to emphasize the different nature of the barrier set in 5 V ). 

Once again one can proceed by separation of variables, finding W(x,t) = 2^^ x 4>n(x)(f>„(A v ) exp(— ^\ n t), with A„ = 
[nn/(A v + Ac)] 2 and <p n (x) = (A v + A c )~ 1/2 sin {\fK.x) . One then gets 

oo / 2 2 \ 

W(M) = A^TX Sm i^+A^) Sm \A^A- X ) ° XP {- 2(A V+ Ay ) ' < A1 > 

n=l N / 

so that one can easily compute the first-crossing rates % and T c across the barriers respectively set in x = A v and x = A c ; 
one obtains 

00 / 2 2 \ 

%{t) - - (XTXfE-- (a7Ta c a ^ cxp {- 2(A n v + Ay ) ^ 

n— 1 



and 



% { t) = J(A V + A c ,t) = TX ^ 2 ^(-D^nsin (^^) cxp [-^fj^A (A3) 

n=l 

respectively. To compute the mass function one only needs /(A, 5 c (z)) = T C (A), which replaced into Eq.(9) gives 

dlnAl^ +1 / hit \ ( n 2 n 2 A(M) 



^ Ad) MM) 



71=1 

However, before going on we need to face a new situation. Both in the case of a single absorbing barrier set at 8 = 8 C 
and in the case of a second reflecting boundary placed at 8 = S v , all the "random-walkers" are eventually going to cross the 
barrier set at 8 — S c ; on the contrary, in the case where the second barrier is also an absorbing one, a relevant fraction of all 
realizations crosses the boundary set at 8 = S v , thereby giving no contribution to collapsed objects. In order to quantify this 
fraction we need to compute the two quantities V v and V c representing the crossing probabilities respectively at 5 = 8 V and 
8 — 8 C . One has 

v.-£ r. W *- f£ (A5, 
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and 

: fV^_ 2^(-l)™ +1 „ : „/ nnA v 
Jo 



J0 n = l 

To compute the series one has to Fourier expand the following elementary functions defined in the interval (0, 2L) 

oo 

... x 11 v-^ 1 . nnx . 
f ( X) =2L = 2-nT,n Sm —> (A7) 



- if0<x < L ^ 
L 2 (-1)™ + . n-KX 

= -> - — sin—, A8 

— - 2 if L < x < 2L n=i 
Lt 

from which one obtains V c = A V /(A V + A c ) and V v = A C /(A V + A c ). One can conclude that, for the cosmologically relevant 
cases, one has V v = S C /(5 V + S c ) = 8 C /(1 + S c ) at any redshift. 

We can interpret this result in two different ways: i) the particles that, because of their peculiar initial conditions, end 
up their existence in "voids" do not form halos at all and then only a fraction V c = 5 V /(S V + S c ) — 1/(1 + S c ) of the total 
mass is collapsed in bound object; ii) the regions that form "voids" are really depleted of mass (this would be obviously 
true if we were following the real dynamics of the field until it reaches the value —1) and then the mass that initially was 
contained in them should have flowed into the overdense regions. The mass function associated with the former interpretation 
is given in Eq. (A4): one of its main properties is that it asymptotically recovers the Press-Schechter form at the high-mass 
end. However, the latter interpretation is the only one that can assure that all the mass is collapsed in object as it is usually 
required in hierarchical models. In fact, if one wants f°° M n(M)dM — {g} to hold, one needs to modify Eq.(A4) and more 
generally Eq. (9) by replacing {q} with the mean density of non-empty regions (qm) = (g) /V c = (1 + 8 c )(g}. 

By introducing the parameter /„ rm (0 < fnorm < S c ) indicating the fraction of mass flowed from underdense into 
overdense regions we can write 



, M , (1 + fnormMQ) A(M) 

n(M,z)= {1+Sc)Hl + z)2 M2 



dlnAlv-^/ ,,„+i . / nn \ ( n 2 -K 2 A(M) \ . 

1 5> 1} nsm (tt^J «p {- ^rrmrhr ) ■ (A9) 

n— 1 x 7 



d\uM\^ y ' \1 + SJ ^\ 2(1 + 5 C ) 2 (1 + ^) 



For scale-free spectra the mass function then becomes 

„ ( „, „ . i^AK (_i£_) - g ( _ ir+ . nsto (_=_) « p (_i£_) . (AI0) 

where M A {z) = 2 1/a M R .(z). 

In Figure Al this solution is compared with the Press-Schechter one: qualitatively the new mass function looks like that 
of Eq.(16). Once again the multiplicity function evolves in a self-similar way; for n p > —2 the maximum value of the mass 
function is reached for M /Ma — [(a/a + l)] 1 ^ 2 ". This implies that this distribution cuts off at a larger mass compared with 
the solution given in the text. 

This feature is present also when one considers the standard CDM power-spectrum (Figure A2); in this case the peak of 
the mass function is reached for M — 5.77 x 10 13 M Q while the typical cluster mass comes out at M c ;(0) = 1.48 x 10 15 M Q . 
The time evolution of our new mass function is consistent with Eq.(18), where (3 ~ 3.27. Once again /3 e ff(z) and a e ff(z) turn 
out to be tightly correlated: their product is approximately equal to one. This agreement is remarkable for z ~ while it gets 
worse as z increases. This is exactly the behaviour we would expect since, as z grows, the integral that defines M c i(z) takes 
contributions from a wider mass interval. 
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Figure Captions 

Figure 1 The mass function n(M) obtained by allowing for mass semi-positivity (solid line) is compared, for different scale- 
invariant spectra, with the Press-Schechter solution (dotted line). The behaviour of the dimensionless and time-independent 
distribution M R n{M)/(g) is shown as a function of M/Mr, where Mr is a redshift-dependent characteristic mass defined in 
the text. 

Figure 2 Time evolution of Mgn(M) / '(g) for a power-law spectrum with n p = 1. 

Figure 3 The present-day mass function obtained by accounting for the mass semi-positivity constraint and the related 
multiplicity function M 2 n(M)/(g) (solid lines) are compared with their Press-Schechter counterparts (dotted lines) for a 
standard CDM scenario. The power-spectrum is obtained starting from a primordial power-law with n p = 1 and using the 
transfer function given by Bardeen et al. (1986) with the choices fix = 1, h = 0.5, as = 1, 5 C = 1.686. 

Figure 4 The mass and the multiplicity functions for the standard CDM model described in Figure 3 are shown at three 
different redshifts. 

Figure Al The time-independent function M\n(M)/(g), obtained by allowing for the existence of void regions and by 
imposing f n0 rm = S c (solid line), is compared, for different scale-invariant spectra, with the Press-Schechter solution (dotted 
line) . 

Figure A2 Present mass and multiplicity functions in a standard CDM scenario. The solutions obtained with the absorbing 
barrier at 5 V and with f n0 rm = 5 C are represented by a solid line, while the Press-Schechter ones are plotted with a dotted 
line. 
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Figure 3 
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Figure 4 
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Figure A2 



